In-Wall Imaging for the Reconstruction of Obstacles by Reverse Time Migration

In this paper, the reverse time migration (RTM) method is applied to the single-frequency reconstruction of embedded obstacles in a wall to perform an introductory study for in-wall imaging. The aim is to determine the geometrical properties of an object embedded in a wall by the use of an information function provided via the RTM method. The method is based on the computation of that information function separately at each point on a reconstruction domain. It is defined as the correlation levels between the incident fields emitted from sources and the back-propagation of the scattered field. The problem is taken from a broader perspective in order to show and confirm the effectiveness of the method. For this purpose, numerical experiments within a fundamental scenario are determined in a particular order to perform an essential Monte Carlo simulation. The paper uses a comparative study to make an objective evaluation of the achievement level of the method in in-wall imaging. The results reveal that the method is at the applicable level of achievement.


Introduction
The reverse time migration (RTM) method based on reverse time extrapolation is proposed in [1]. Besides many other migration algorithms [2], it is proven to be of high accuracy in the reconstruction of complex structures [3][4][5]. Among direct methods including [6][7][8][9][10][11] and many others, the RTM method was also used and useful for imaging obstacles in several applications. Studies were conducted using acoustic waves within a free space application in [12], within a planar waveguide in [13], within a half-space in [14], with only intensity measurement in [15]; using electromagnetic waves within free space in [16], within a rectangular waveguide in [17], within an application based only on the intensity of data in [18], within a biomedical tomography at optical frequencies in [19]; and using elastic waves within a half-space in [20].
Applications governed by the scalar wave equation are more feasible due to computational expenses arising as a challenging problem. Employing an asymptotic filter to map the real 3D data, which are obtained on an acquisition surface far enough from the reconstruction domain into a 2D space, creates extra room to reduce that cost to lower levels [23,24]. The spatial conversion filter approach in [23] also underlines the waveform, i.e., the argument of a received signal, whereas the phaseless magnitude approach is of concern in [15,18]. Consequently, each part of the data is proven to have its own benefit.
Imaging of sparse objects with the Born approximation and the compressive sensing approach is studied in [36]. It was proposed that the hypothesis of sparsity makes it possible to reduce the size of the data. A novel approach using the valuable features of the linear sampling method is presented in [37] for imaging the geometrical properties of an embedded object. A microwave tomographic approach is employed to process ground penetrating radar (GPR) data for detecting and locating defects in a historical wall [38]. A shadow projection method is studied to determine the geometrical properties of air gaps in a dense medium [39]. The shadow projection method is based on the data obtained on the boundary behind the object. The method exploits the lensing effect emerging between the two adjacent layers with high contrast. A prototype based on impulse radio ultra-wideband radar is developed for identifying the current condition of a wall [40]. The prototype includes a deep learning module for a robust and self-adaptive application. Synthetic aperture radar imaging is employed for the reconstruction of the infrastructure elements within a dry wall in [41]. For this purpose, a frequency-modulated continuous-wave radar operating at 80 GHz with a bandwidth of 25.6 GHz is used.
The most basic wall parameters are dielectric permittivity, conductivity, and thickness. Reconstruction of wall parameters is conducted for several reasons such as determining the thickness of a wall, checking its health condition, detecting the infrastructures within, and planning an indoor wireless broadcast scheme [42][43][44][45][46][47].
The signal response of a wall is always stronger than that of a passive object located behind it and therefore a clutter reduction procedure is needed in TWI applications [48][49][50][51].
Based on the RTM method, a radar application is conducted for imaging targets buried in a multilayered underground with the help of probes located inside the test site through a dug hole [32]. For the sake of the structural reliability in some components, imaging of notches and defects on aluminum samples by using an ultrasonic array leaning on the material at zero distance to excite elastic waves is introduced in [33,34]. The internal defect detection of casting steel is studied in [35].
Nondestructive testing techniques constitute a big family. One member of this family is microwave imaging, which uses electromagnetic waves in the microwave frequency band [52][53][54]. Electromagnetic illumination is achieved by using this frequency band throughout this study.
The main aim of this paper is to produce fundamental results for in-wall imaging with the RTM method for the detection and imaging of embedded obstacles at microwave frequencies. Within this context, imaging an embedded object in the inaccessible middle layer of a three-layered medium is considered.
The rest of the study presented here is organized as follows. In Section 2, the geometry of the problem and the general formulation are given. In Section 3, the RTM method based on obtaining support within a reconstruction domain is presented and the application guidelines of the Monte Carlo simulation within this study are introduced. In Section 4, a comprehensive application through a properly defined set of numerical experiments is conducted to observe the performance of the algorithm for in-wall imaging. Additionally, the outline of the quantitative performance evaluation is included and discussed in this section. Section 5 is the Conclusion section revealing a summary as concluding remarks for the study.
Wave excitation is achieved under the time-harmonic regime. To this end, a time factor e −iωt is selected to be dropped in this paper including the whole numerical experiments conducted and presented in Section 4.

Formulation of the Problem
The reference geometry of the problem considered here is a wall structure of finite thickness laying along the horizontal axis, basically a three-layered medium configuration with the object of interest sitting on a finite portion of it. That geometry is depicted in Figure 1. The wall is defined as a lossless material with a relative permittivity ε r 2 = ε r(wall) and a thickness |d 1 − d 2 |. The first and the third layers are assumed to be free space, ε r 1 = ε r 3 = 1. The magnetic permeability is the same as that in free space, µ r = 1, in every area. The reconstruction domain, denoted by Ω, is a subset of the wall being the region of interest that the object (domain), D, is located within. Time-harmonic field excitations are achieved by infinite line sources and a 2D geometry is considered. Locations of the source and the observation points are selected to be identical and they are positioned on both or one of the two sides of the wall forming a multibistatic configuration for data acquisition, which is denoted as Γ.

Formulation of the Problem
The reference geometry of the problem considered here is a wall structure of finite thickness laying along the horizontal axis, basically a three-layered medium configuration with the object of interest sitting on a finite portion of it. That geometry is depicted in Figure 1. The wall is defined as a lossless material with a relative permittivity = ( ) and a thickness | − |. The first and the third layers are assumed to be free space, = = 1. The magnetic permeability is the same as that in free space, μ = 1, in every area.
The reconstruction domain, denoted by Ω, is a subset of the wall being the region of interest that the object (domain), , is located within. Time-harmonic field excitations are achieved by infinite line sources and a 2D geometry is considered. Locations of the source and the observation points are selected to be identical and they are positioned on both or one of the two sides of the wall forming a multibistatic configuration for data acquisition, which is denoted as Γ. Wave dynamics in the presence of an infinite line (a TM polarized) source is described by the Helmholtz wave equation as: where = ( , ) denotes the electric field deviation; is the wavenumber; is the operating frequency; is the magnetic permeability in the medium; is the amplitude of the electric current source; ( , ) is the cartesian coordinate system denoting points of field (observation); and ( , ) is the point of a line source excitation [55]. The Helmholtz equation is described as a differential operator acting on a field distribution and the solution is given by the (scalar) Green's function, , as: resulting in a fundamental solution to Equation (1) under the conditions determined due to the boundaries of the layers at a problem geometry of interest as: and are continuous on the boundaries, with the radiation condition as the distance from the excitation point, | ⃗ − ⃗ |, tending to infinity [55]. Here, ⃗ = ( , ) ≡ ( , ) denotes the polar coordinate system and ⃗ = ( , ) ≡ ( , ) is the point of the line source excitation.
The solution of Equations (2) and (3) for ( , ; , ) is given as an explicit expression due to the problem geometry and is explained in detail in the literature for many cases including free space and multilayered media geometries [55][56][57][58]. Under the time convention given in this paper, the Hankel function of the first kind, ( ) , is reserved for Wave dynamics in the presence of an infinite line (a TM polarized) source is described by the Helmholtz wave equation as: where E = E(x, y) denotes the electric field deviation; k is the wavenumber; ω is the operating frequency; µ is the magnetic permeability in the medium; I is the amplitude of the electric current source; (x, y) is the cartesian coordinate system denoting points of field (observation); and (x 0 , y 0 ) is the point of a line source excitation [55]. The Helmholtz equation is described as a differential operator acting on a field distribution and the solution is given by the (scalar) Green's function, G, as: resulting in a fundamental solution to Equation (1) under the conditions determined due to the boundaries of the layers at a problem geometry of interest as: G and ∂G ∂y are continuous on the boundaries, with the radiation condition as the distance from the excitation point,  (2) and (3) for G(x, y; x , y ) is given as an explicit expression due to the problem geometry and is explained in detail in the literature for many cases including free space and multilayered media geometries [55 -58]. Under the time convention (2) 0 , is for incoming ones. In free space, Green's function expression is given as follows: where k 0 denotes the wavenumber of the free space.
In a three-layered medium, the problem geometry is more complicated and the expression of G(x, y; x , y ) occupies considerable space. To this end, all expressions are given in Appendix A.
The pair of scattering equations for the field scattered from a dielectric object, E s , due to an incident field, E i , are described by those two Fredholm integral equations as [55][56][57][58].
where (x, y) are the points on a selected acquisition line; (x , y ) are the points on the object domain or the reconstruction domain; and χ =ε r (x, y) − 1 is the object's function withε r denoting the complex permittivity. k is the wavenumber of that layer at which the object domain is located. E is the total electric field inside the reconstruction domain and is described by the implicit expression in Equation (6) [55 -58]. Similarly, the field scattered from a perfect electric conductor (PEC) object is described by another pair as follows [55-58]: Here, (x , y ) are the points located on the surface of the PEC object. ∂D denotes the boundary of the object domain, D. j and E denote the inducted electric current deviation and the total electric field on that surface, respectively. The expression given in Equation (8) is, in fact, another form of Equation (3) due to the presence of a PEC discontinuity [55-58].

Reverse Time Migration (RTM) Method
By obtaining the back-propagation and computing the cross-correlation, the RTM method employed here is based on creating support for obstacles embedded in a medium by using the far-field scattering data, E s , [12][13][14][15][16][17][18][19][20]. It is assumed that there are N s and N r locations on Γ s = ∂B s and Γ r = ∂B r as points of source and observation, respectively. B s and B r denote continuous line segments, open or closed curves where N s points of source and N r points of observation are located at some points of them, respectively: Γ s,r = ∂B s,r ⊂ B s,r . The reconstruction domain may fall inside both or one of those two with obstacles in it, or may not.
The back-propagation, E b , is equivalent to the propagation excited by the timereversed form of the scattered data which is the complex conjugation in the frequency domain. The time-reversed scattered data stand as a source term in the Helmholtz equation that will be solved for back-propagation [12][13][14][15][16][17][18][19][20]: Sensors 2023, 23, 4456 Here, (x s , y s ) and (x r , y r ) denote the points of source and observation on Γ s and Γ r , respectively, and (x, y) denotes any point on a continuous reconstruction domain. |Γ r | N r is placed in Equations (9) and (11) for stabilization and |Γ r | denotes the length of Γ r . u s and u b are the fundamental expressions of the scattered field and the back-propagation, respectively: E s = iωµI·u s and E b = iωµI·u b . u s denotes the complex conjugation of u s .

Monte Carlo Simulation
Reporting the level of achievement of the RTM method in in-wall imaging based on a single test would not be fair and would also be misleading. In order to make a general inference, the method should be tested under different circumstances. The Monte Carlo simulation scheme can be adapted for this purpose [59]. A properly defined set of experiments may be employed in the outline of this approach within a comparative study. Each experiment results in a different level of achievement and an overall evaluation will be possible by combining all of them together. Then, it is possible to make a proper conclusion.
In order to achieve this goal, first, a set of parameters is defined. Those parameters represent the domain of possible inputs in the Monte Carlo simulation scheme and they are listed as follows:
Material type of the embedded object; 3.
Structural property of the embedded object; 4.
Location of the embedded object; 5.
Noise level in the medium.
Any change in any parameter results in a different experiment with different scattered data. The scattered data, u s (x r , y r ; x s , y s ), is synthetically obtained from the pair of scattering equations given in Equations (5) and (6) for a dielectric object and in Equations (7) and (8) for a PEC object. Hence, the noise level in the medium is added to the scattered data. The Green's functions are G(x, y; x s , y s ) and G(x, y; x r , y r ) and they are synthetically obtained from Equation (4) or (A1)-(A13). In the Monte Carlo simulation scheme, the scattered data and Green's functions are taken instead of randomly generated inputs. The deterministic computation step is defined in Equation (12). After all of those, the level of achievement is calculated quantitatively for each experiment to make a final assessment.
The block diagram of the Monte Carlo simulation scheme adapted for the RTM method is given in Figure 2. cally obtained from Equation (4) or (A1)-(A13). In the Monte Carlo simulation scheme, the scattered data and Green's functions are taken instead of randomly generated inputs. The deterministic computation step is defined in Equation (12). After all of those, the level of achievement is calculated quantitatively for each experiment to make a final assessment.
The block diagram of the Monte Carlo simulation scheme adapted for the RTM method is given in Figure 2. The blocks of the diagram in Figure 2 are described as follows: • Input parameters are the six parameters listed above; • I ori is the original information data; • G(x, y; x s , y s ) and G(x, y; x r , y r ) are the Green's functions computed on (Ω; Γ s ) and (Ω; Γ r ), respectively; • If the embedded object is a dielectric material, Green's functions are computed on two different domains, G(x, y; x , y ) on (Γ r ; D) and (D; D); • If it is a PEC material, the inducted electric current is computed on the surface of the object, j(x , y ) on ∂D; • u s (x r , y r ; x s , y s ) is the scattered data synthetically acquired on Γ r due to the sources located at Γ s ; • The noise level in the medium is added to the noise-free data: u s (x r , y r ; x s , y s ) + N; • I rec is the reconstructed information data obtained from the RTM method; • Procedures for quantitative performance evaluation of an experiment and an overall assessment of the whole experiment set are considered in Section 4.

Numerical Experiments and Discussion
In order to show the performance of the algorithm in line with reconstructing objects embedded in a wall, part of which is equivalent to those applications given in [12,16]. In other words, those examples of adaptation to in-wall imaging are considered by copying some of the parameters except those determined by the problem geometry. For instance, points of source and observation cannot be located inside the wall within the sense of complete nondestructive testing and using an additional probe inside the wall is not considered. Thus, it will be possible to establish a comparative study.
Covering the range y = (−3, +3) m, there is a wall of thickness at 6 m and it lies along the whole Ox axis. The relative permittivity of the wall is ε r(wall) = 2. The reconstruction domain covers some limited part of that wall as it extends within the range of Ω = (−3, +3) m × (−3, +3) m. The domain, Ω, is laid out on a grid pattern of 200 × 200 meshes of equal size and the correlation functionals are computed at their central points [56].
A set of experiments, each touching upon a different aspect of the subject, is created to acquire different performance outputs of the method due to different situations.
The first parameter is the acquisition line. Four different sets of probes are used: In the second and third sets, Γ 2 and Γ 3 , are formed by using only the upper half and the lower half of the first set located at the same line segments at y = +4 m and y = −4 m, respectively: Γ 2,3 = ∂B 2,3 and N 2,3 = 97. One-sided illumination and measurement is performed in case one of the two sides of a wall is inaccessible; • There is also a zeroth set, Γ 0 = ∂B 0 , formed by different points of source and observation, which are uniformly distributed on a circle with radius r = 12×8 π ≈ 5.53 m and N 0 = 194 illuminating the same size of the area as that in the first set and it is given only in the absence of the wall for comparison.
The reconstruction domain, Ω, is surrounded by both the first and the zeroth sets and sits in the middle of both those two different closed curves, B 0 and B 1 .
Points of fully isotropic source and observation uniformly located on a circle with a radius big enough to encapsulate the reconstruction domain in an unbounded medium from all directions (here, it is Γ 0 ) stands as a reference level of perfection for data acquisition and using another acquisition line instead mostly leads to some imperfections in reconstruction [16]. In a non-destructive in-wall imaging application, there are some restrictions on determining such a proper acquisition line as it lacks that kind of perfect surrounding property. There cannot be any penetration into the wall, therefore it ends up with some probable imperfection, as expected [16].
The second parameter is the material type. The object may be selected from three different types of materials: • First, a dielectric object of penetrable material with a contrast ratio of ε r(object) /ε r(medium) = 2; • Additionally, an air gap only in the presence of the wall, again, as a penetrable material, with a contrast ratio of ε r(object) /ε r(wall) = 0.5; • Thirdly, a PEC object of non-penetrable material.
The third parameter is the structural property of the embedded object. Based on this, two canonic geometries are selected:

•
A unit circle, bigger in size carrying a basic curvature property; • A point body that stands for smaller objects being poor in structural properties.
The fourth parameter is the location. The object may be located at four different positions: The sixth parameter is the noise level and four different noise levels are selected: N ∼ 0%, 10%, 20%, and 50%.
There were 1632 numerical experiments conducted. In total, 768 of them were free space applications and the remaining 864 were conducted in the presence of a wall. There were 96 more than the initial 768 since it is also another option to locate or embed material into a wall producing contrast in the reverse order, i.e., a material of some permittivity less than that of the medium, which here was an air gap. One could find out that each of the 192 experiments of those 768 is conducted with Γ 0 , Γ 1 , Γ 2 , and Γ 3 , and each of those 288 from the remaining 864, which is again 96 more, is carried out with Γ 1 , Γ 2 , and Γ 3 , respectively. Any change in any parameter reveals a different experiment. By combining them all together, an essential multiparameter Monte Carlo simulation is established through all those experiments for the assessment of an overall fundamental performance evaluation of the method in in-wall imaging applications.
The results for some of those cases in which the objects are located at the most unbalanced position are given in Figures 3 and 4 for the dielectric and the PEC unit circles, and in Figures 5 and 6 for the dielectric and the PEC point bodies, respectively. The plots are given only for the imaginary parts of the correlation functionals [12][13][14][15][16][17][18]20] due to sign compliance [12,13,16]. In Figure 3, in the first row, the dielectric unit circle reconstruction with Γ 0 is given in the absence of the wall for comparison. Reconstruction of the same object with the double-sided multibistatic configuration, Γ 1 , again, in the absence of the wall for comparison and also in the presence of it are given in the second and the third rows, respectively. Reconstruction of the air gap is given in the fourth row. The PEC unit circle reconstructions are given in Figure 4 in the same order as in Figure 3. The point body reconstructions are given in Figures 5 and 6 in the same order as in Figures 3 and 4.
Such contour plots placed side-by-side, as those given in Figures 3-6, present instances for visual comparisons [60]. However, it is hard to state that visual comparison is suitable for an overall evaluation of the whole of the 1632 experiments and the outputs. Mean squared error (MSE) of the support functions, I rec , with respect to the original data, I ori , may facilitate an efficient quantitative evaluation. Therefore, a proper definition must be made for modeling the original data. A mesh on a given grid pattern with more than 50% of it covered by part of the object is valued at 1, otherwise, it is 0. Additionally, a PEC object does not allow electromagnetic waves to penetrate inside the region bounded by itself and therefore they are stuffed materials whereas dielectric objects are penetrable. Therefore both PEC circles and disks of the same size are taken as stuffed and identical to each other. On the other hand, dielectric circles and disk structures are different due to wave penetration and a dielectric circle is taken as empty or hollow. Suffering from dimension, there is no such discussion on point bodies.
Working with raw data for the calculation of MSE values is similar to comparing apples with oranges. Using the min-max scaling, a normalization procedure is followed for linearly mapping the values of the support functions into the range [0, 1]. Denoting the normalized data asÎ rec , the MSE calculation for I ori andÎ rec instead of I ori and I rec is carried out and the quantitative performance evaluation percent, P, is given in the outline of this definition, accordingly: P(I ori , I rec ) = 1 − MSE I ori ,Î rec × 100. with Γ is given in the absence of the wall for comparison. Reconstruction of the same object with the double-sided multibistatic configuration, Γ , again, in the absence of the wall for comparison and also in the presence of it are given in the second and the third rows, respectively. Reconstruction of the air gap is given in the fourth row. The PEC unit circle reconstructions are given in Figure 4 in the same order as in Figure 3. The point body reconstructions are given in Figures 5 and 6 in the same order as in Figures 3 and 4.    It is important to note that as the area of a reconstruction domain tends to infinity, any embedded object becomes a truly dimensionless point body and the MSE value approaches a perfect value of 0. This means that the MSE calculation can only be of use in a comparative study within a finite reconstruction domain in order to make a proper evaluation. As a consequence of this, it may be expected that the MSE calculation would be in favor of the point body reconstructions instead of those conducted with a unit circle since a unit circle is bigger in size than a point body. On the other hand, bigger objects return greater responses on the same acquisition line. The magnitude values of a support function on those meshes at which part of the object is located are mostly much greater than the rest and one may look that situation up in Figures 3-6. As a result of those two, the formulation proposed in (13) becomes reasonable for a quantitative performance evaluation.
In line with the parameters listed above, the Monte Carlo simulation outputs are given to show the performance levels by acquisition line, material property, structural property, object location, operating wavelength, and noise level. The results are given in a crosstab form in Tables 1-6, respectively. In the first row of the tables, those experiments conducted with the acquisition line Γ 0 are given in the absence of the wall for comparison. The reconstructions with Γ 1 , Γ 2 , and Γ 3 , again, in the absence of the wall for comparison and also in the presence of it are given in the second and the third rows, respectively. On the far right, the overall evaluation is always placed in bold issuing a condensed brief of the study.
(i) (j) (k) (l) Figure 4. Reconstructions of the PEC circle are given in (a-l) in the same order as in Figure 3.  In Table 6, the robustness of the algorithm is shown due to different levels of additive uniform noise, which is also presented for visual comparison in  (m) (n) (o) (p) Figure 5. Reconstructions of the dielectric and the air gap point bodies are given in (a-p) within the same order as in Figure 3. Such contour plots placed side-by-side, as those given in Figures 3-6, present instances for visual comparisons [60]. However, it is hard to state that visual comparison is suitable for an overall evaluation of the whole of the 1632 experiments and the outputs. Mean squared error (MSE) of the support functions, , with respect to the original data, , may facilitate an efficient quantitative evaluation. Therefore, a proper definition must be made for modeling the original data. A mesh on a given grid pattern with more than 50% of it covered by part of the object is valued at 1, otherwise, it is 0. Additionally, a PEC object does not allow electromagnetic waves to penetrate inside the region bounded by itself and therefore they are stuffed materials whereas dielectric objects are penetrable. Therefore both PEC circles and disks of the same size are taken as stuffed and identical to  The presence of a wall has veiling and bending effects on electromagnetic waves. As an electromagnetic wave hits a wall, some portion of it does not penetrate into it to reflect back, resulting in a reduction in signal power in the first place. On the other hand, the wall triggers multiple reflections and transmissions in the medium and a series of interactions emerges between the (boundaries of the) wall and the embedded object increasingly adding some extra power to the signal strength. This situation is also seen in the value scales of the contour plots in Figures 3-6. Overall results show that, in most cases, those interactions gain an advantage and have constructive effects on reconstruction performance in the end.

Conclusions
A correlation-based algorithm, the RTM method, is studied for the reconstruction of objects within an inaccessible wall. A comparative study is presented. Numerical experiments put forth essential achievements giving promising results for further studies and making it feasible for in-wall imaging applications.
Inevitably, using a set of probes located outside of the wall instead of perfectly surrounding the reconstruction domain deforms the reconstruction. Interactions between the embedded obstacle and the boundaries due to the presence of a wall geometry occur to work as a retouching tool in the horizontal direction and have positive effects on the reconstruction. On the other hand, those interactions produce interferences deforming the reconstruction in the vertical direction. As the distance of an object away from the boundaries of the wall increases/decreases, it also improves/deteriorates the imaging quality. With a decrease in operating wavelength or at higher operating frequencies, that distance is relatively on the rise.
Much more complicated problem geometries and real-life applications may be taken into consideration for future works in order to see where the efficiency and the applicability levels of the method could reach further in in-wall imaging.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: RTM reverse time migration TWI through-wall imaging 3D three-dimensional 2D two-dimensional GPR ground penetrating radar TM transverse magnetic PEC perfect electric conductor MSE mean squared error SIP Sommerfeld integration path